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We present the first numerical calculation of the (local) metric perturbation produced by a small 
compact object moving on an eccentric equatorial geodesic around a Kerr black hole, accurate to 
first order in the mass ratio. The procedure starts by first solving the Teukolsky equation to obtain 
the Weyl scalar ip 4 using semi-analytical methods. The metric perturbation is then reconstructed 
from ip 4 in an (outgoing) radiation gauge, adding the appropriate non-radiative contributions arising 
from the shifts in mass and angular momentum of the spacetime. 

As a demonstration we calculate the generalized redshift U as a function of the orbital frequencies 
fl r and fl# to linear order in the mass ratio, a gauge invariant measure of the conservative corrections 
to the orbit due to self-interactions. In Schwarzschild, the results surpass the existing result in the 
literature in accuracy, and we find new estimates for some of the unknown 4PN and 5PN terms in 
the post-Newtonian expansion of U. In Kerr, we provide completely novel values of U for eccentric 
equatorial orbits. Calculation of the full self-force will appear in a forthcoming paper. 


I. INTRODUCTION 

With Advanced LIGO scheduled to start data collection 
this year T], the detection of gravitational waves seems 
imminent. One of the main expected sources of gravi¬ 
tational waves are binary systems of compact objects in 
close orbits. Detection of gravitational waves from these 
systems is dependent on our ability to model them ac¬ 
curately, requiring us to solve the two-body problem in¬ 
cluding general relativistic effects. Unlike its Newtonian 
counterpart, the general relativistic two-body problem 
admits no analytic solutions and requires some sort of ap¬ 
proximation scheme to solve. There are currently three 
main approaches, which are largely complementary. 

Numerical relativity (NR) works by discretizing space- 
time and numerically solving the nonlinear field equa¬ 
tions of general relativity on a grid. Spectacular progress 
has been made in this approach over the last few decades 
(see [2] and the references therein). However, discretiza¬ 
tion of spacetime means that this approach has lim¬ 
ited use in systems with vastly differing length and time 
scales. 

Post-Newtonian (PN) theory (effectively) expands the 
equations of motion in the binary separation, obtain¬ 
ing corrections to Newtonian theory order-by-order (see 
[3j and the references therein). This provides extremely 
good models for the motion at large separations. But as 
one enters the strong gravity regime, convergence of the 
post-Newtonian expansion is poor. 

Black hole perturbation theory (used in this paper) ap¬ 
proximates the motion of the binary by expanding in the 
mass-ratio of the components, starting from test par¬ 
ticle (geodesic) motion in the geometry generated by 
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the largest mass[3HB]. This is most naturally applied 
to systems with a very small mass ratio, such as ex¬ 
treme mass ratio inspirals (EMRIs) formed by a stellar 
mass compact object spiralling towards a (super)massive 
black hole, expected to occur regularly in galactic cores. 
The gravitational waves produced by an EMRI fall out¬ 
side the sensitivity range of ground based interferometers 
(LIGO/Virgo/KAGRA), but are expected to be one of 
the prime sources for the space-based gravitational wave 
observatory eLISA, slated for launch in the 2030s [7]. 

In the meantime, black hole perturbation theory can be 
used to validate and inform NR and PN techniques by 
comparing results in regions of the compact binary pa¬ 
rameter space where its range of applicability overlaps 
with the other methods. To compare results between 
different approximation methods requires studying quan¬ 
tities that are insensitive to the approximation used. In 
particular, since different approaches tend to work in dif¬ 
ferent gauges, we need the quantities to be gauge invari¬ 
ant. 

In addition, results can be used to tune the effective-one- 
body (EOB) formalism developed by [8]. This formal¬ 
ism takes results from PN theory, black hole perturba¬ 
tion theory and numerical relativity to reduce the two- 
body problem to an effective analytical model of a single 
point particle moving in an effective metric j9j similar to 
that of the Newtonian reduced mass problem. Advances 
have been made with the latest being the evolution of 
a non-spinning binary-neutron-star by mm whose re¬ 
sults agree to those of NR within numerical accuracy. 
The free parameters in this theory are fixed by compar¬ 
ing gauge-invariants of physical interest from the above 
three approaches. 

First such invariant to be calculated was the “redshift” 
invariant |T7| , the linear-in-mass-ratio change in the time- 
component of the 4-velocity of the compact object in cir¬ 
cular orbit around a Schwarzschild black hole which was 



2 


then used to calibrate the binding energy and angular 
momentum of a binary in the EOB model Il3l . Later 
on, many other invariants have been calculated for cir¬ 
cular orbits, including corrections to the circular limit of 
the periastron advance j!4j. change in innermost-stable- 
circular orbit p51 \W \, spin precession E, quadrupolar- 
tidal [l 8 j , and octupolar-tidal invariants m ■ PN expan¬ 
sions of the latter three have been studied to very high 
order IBS ED], and were used to calibrate the EOB model 
further [ 21 ; . 

Calculations of the redsliift invariant for circular orbit 
in Schwarzschild spacetime have been generalized to ec¬ 
centric orbits around a Schwarzschild black hole where 
they calculated the orbit-averaged value of the invari¬ 
ant, U, for given azimuthal and radial frequencies using a 
Lorenz gauge in the time-domain |22| , and recently, a suc¬ 
cessful comparison was performed for this orbit-averaged 
quantity between the self-force and PN theories using 
frequency domain methods [23] . 

All these calculations for orbits around a non-rotating 
(Schwarzschild) black hole have been performed either 
in a Lorenz gauge or in the Regge-Wheeler-Zerilli gauge 
(or both). In both gauges, the linearized Einstein equa¬ 
tion is separable. In the frequency-domain, they may be 
solved mode-by-mode by solving second order ordinary 
differential equations. 

A major obstacle in generalizing these calculations to 
rotating (Kerr) black holes has been that there are no 
known gauges for which the linearized Einstein equation 
in Kerr spacetime is separable. In recent years, an alter¬ 
native approach has been put forward that exploits the 
fact that the Teukolsky equation for the Weyl scalars 4>o 
and ^4 are separable and can be solved efficiently. Using 
the formalism developed by Chrzanowski [21] , Cohen and 
Kegeles [251126], and Wald [271 (CCK or CCKW), one can 
reconstruct the metric perturbation in a radiation gauge 
from the Weyl scalars. 

This approach has been implemented to obtain the met¬ 
ric perturbation produced by a particle moving on a 
circular orbit around Schwarzschild [281 25] , and Kerr 
black holes [3D]. Using this technique high precision re¬ 
sults for some of the invariants above have been obtained 
in the Schwarzschild case and employed to extract very 
high PN-expansion using a numerical fitting technique 
unusmm In the Kerr case, the shift in the inner¬ 
most-stable circular equatorial orbit fj], and the red- 
shift invariant (30] have been calculated. Work on the 
PN-expansion of the latter is in progress [53] . 

Despite the success of [25H3CT] , it was unclear that the self¬ 
force corrections to the motion that they calculated were 
well-defined due to radiation gauge metric perturbations 
being irregular in a neighbourhood of the particle world¬ 
line. The use of radiation gauge metric perturbations to 
obtain self-force corrections to the motion of a particle 
was put on a firm formal footing in [34] . In particular, 


they showed how the mode-sum regularisation formula 
for the self-force needs to be modified to use radiation 
gauge metric perturbations as its input. 

The goal of this paper is to implement the formalism set 
out in [33] for eccentric (equatorial) orbits to obtain— 
for the first time—the metric perturbation produced by 
a particle moving on an eccentric orbit around a Kerr 
black hole. A problem with extending this approach 
to eccentric orbits is that the CCK procedure for ob¬ 
taining the metric perturbation from the Weyl scalars is 
only well-defined for vacuum spacetimes. However, when 
solving the Teukolsky equation sourced by a particle on 
an eccentric orbit in the frequency domain, the particle 
point source gets smeared out of the entire region be¬ 
tween the periapsis r min and apapsis r max of the orbit. 
We overcome this problem utilizing the method of ex¬ 
tended homogeneous solutions [35] to avoid dealing with 
non-vacuum solutions at any stage of the calculation. 

Our procedure will be as follows. We obtain highly accu¬ 
rate solutions to the Teukolsky equation using the semi- 
analytical Mano-Suzuki-Takasugi (MST) formalism. We 
then algebraically invert the fourth order differential 
equation for the intermediate Hertz potential, from which 
we obtain the (singular) metric perturbation using the 
CCK procedure. The regular part of the metric pertur¬ 
bation is obtained using mode-sum regularization. The 
missing pieces of the metric perturbation due to pertur¬ 
bations of the mass and angular momentum of the back¬ 
ground spacetime, which cannot be recovered using the 
CCK procedure, are obtained using the result from (SB] , 

Using the obtained (regular part of the) metric pertur¬ 
bation we calculate the self-force correction to general¬ 
ized redshift invariant U for eccentric orbits as a function 
of the orbital frequencies. In Schwarzschild, our results 
match to the previously published results of [22] and [23] 
to all given digits, surpassing them in accuracy and com¬ 
putational efficiency. By fitting to a large dataset of or¬ 
bits, we are able to recover all known coefficients of the 
PN expansion of the generalized redshift in Schwarzschild 
and obtain estimates for some of the unknown 4PN and 
5PN terms. 

In Kerr, our result matches the previously published in 
circular limit [30]. Our completely novel results for the 
generalized redshift for eccentric equatorial orbits around 
a Kerr black hole, pass all the consistency checks that we 
can perform. In particular, we numerically match the an¬ 
alytically known leading regularization parameters. Be¬ 
sides verifying our numerical method, this provides the 
first every test of their analytical calculation for eccentric 
orbits in Kerr. 

The paper is organized as follows. In section [TT] we re¬ 
view some of the background need to setup our calcu¬ 
lation, introducing the notation and conventions we will 
need. In particular, we will define the generalized red¬ 
shift invariant that we will be calculating. Section m 
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then describes the method that we will use to calculate 
the metric perturbations generated by a particle orbiting 
a Kerr black hole on an equatorial eccentric orbit. Details 
of the numerical implementation of this method are given 
in section [TV| In section [V] we numerically calculate the 
generalized redshift invariant U. We first compare the 
results from our code to values in the existing literature 
in limits of eccentric orbits around a Schwarzschild black 
hole and circular equatorial orbits around a Kerr black 
hole. We then give the first ever numerical values of U for 
eccentric equatorial orbits around a spinning black hole. 
We conclude in section fVTl with a discussion of the limita¬ 
tions of our code, and how it can be extended to calculate 
different quantities derived from the metric perturbations 
including the full (first order) self-force correction to the 
motion. 


II. PRELIMINARIES 


A. Some conventions and notation 


In this paper, we use geometric units c = G = 1. In 
addition, we will usually set the mass of the central black 
hole, M to 1 as well. 

Whenever a sum is written without explicit bounds on 
the indexed summed over, it is assumed to be summed 
over its full natural range, i.e. all (integer) values of the 
index for which the summand is well-defined and non¬ 
zero. For example, sums over £ will typically range from 
\a\ to oo, sums over m will range from — £ to £, etc. 

When in this paper we refer to (spin-weighted) spheri¬ 
cal/spheroidal harmonics, we generally mean the “polar” 
part of the function, i.e. the part that depends on the 
polar coordinate 9 (or z , see below). These functions are 
normalized such that 



d z s Y lm {z) 2 = 1. 


The usual spherical harmonics are therefore 


Ylm(OA) 


O y im (cosfl)e*"* 

V2tt 


( 1 ) 

( 2 ) 


We use overbars to denote complex conjugation, i.e. x is 
the complex conjugate of x. 

In this paper, we write the metric generated by a rotating 
black hole with mass M = 1 and spin a as, 


2 r 


= - (! - -£-) d * 2 + ^ dr2 + r 


dz 2 


1 -z 


+ — /^(2a 2 r(l-z 2 ) + (a 2 +r 2 )£)d</> 2 (3) 


4ar(l — z 2 ) 


dt d <f>, 


with 


A = r(r — 2) + a 2 , (4) 

£ = r 2 + a 2 z 2 . (5) 

Here t, r, and (f> are the usual Boyer-Lindquist coordi¬ 
nates, while z is related to the usual Boyer-Lindquist co¬ 
ordinate 9 by z = cos 9. In particular, the equator of our 
black hole spacetime is given by z = 0. The two roots of 
A indicating the location of the inner and outer horizon 
are denoted 


r± = 1 ± \/l — a 2 . (6) 


Much of the techniques used in this paper rely on the 
Newman-Penrose formalism. We use the following null 
tetrad (the Kinnersley tetrad in our modified Boyer- 
Lindquist coordinates), 


^ = 1“ = -^(r 2 + a 2 , A,0,a), 


= «T = —(r 2 +a 2 ,-A,0,a), 


eii = m M = — 1 


p\J 1 — z 2 


V 2 


e(f = = 


py/l — z 1 


V2 i_ z2 )> 


with 


-1 


P = 


(7) 

( 8 ) 

(9) 

( 10 ) 


( 11 ) 


We will use Greek letters to represent spacetime indices 
and Latin letters for tetrad indices. In particular, 


gab — £a9^v e b 


/ 0 -1 0 0 \ 

I -1 0 0 0 1 
0 0 0 1 ’ 
Vo 0 10 / 


( 12 ) 


is the metric in tetrad indices, and will be used to raise 
and lower tetrad indices. 


B. Two-body problem in black hole perturbation 
theory 


We are interested in the motion of a pair of gravitation¬ 
ally bound compact masses m and M, in the limit that 
m -C In this limit, the gravitational field produced 
by the smaller object m can be treated as a perturbation 
to the black hole geometry generated by the larger object 


1 With some abuse of notation we will use m and M both as labels 
for the two objects and as the value of their masses. 
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M. If we place ourselves in a frame in which M is at rest, 
we can study the motion of m as a perturbative series in 
the mass ratio p = m/M <C 1. 

Both masses are assumed to be compact enough that 
they can be treated as black holes. Moreover we assume 
both masses to have no charge. We allow the larger mass 
to have non-zero spin a = J/M. The background geom¬ 
etry is therefore described by the Kerr metric, which in 
addition is assumed to be sub-extremal, i.e. a/M < 1. 
We further assume the smaller object to have zero spin. 


C. Geodesics in Kerr spacetime 


In the limit p —> 0, the smaller object will follow a 
geodesic of the Kerr spacetime generated by the larger 
object M. The geodesic equations in Kerr spacetime have 
the following form: 


^dr^ 2 

R(r) 

(13) 

Vdr/ 

E(r,z) 2 ’ 

^dz^ 2 

Z(z) 

(14) 

Vdr/ 

S(A^) 2 ’ 

d (j) 

dr 

$ r (r) + $ 2 (z) 
E(r,z) 

(15) 

d t 

dr 

T r (r) + T z (z) 

E(r,z) 

(16) 


where R, Z , <b r , <f> z , T r , and T, are known functions of 
r and z (see e.g. m)- 

As first shown by Carter (38], this set of equations is 
separable. This separation can be achieved easily, by 
changing to a convenient time variable A to parametrize 
the orbit, 


S^s=E(r,*) =r 2 + a 2 z 2 . (17) 

dA 

The time variable A is commonly referred to as “Mino 
time”. This parametrization allows for an analytic solu¬ 
tion of geodesic equations in terms of elliptic functions 

EH- 

The solutions in |37) are most naturally expressed in 
terms of the periapsis r m - ln , the apapsis r max , and the 
maximum value of the z coordinate z max . In this paper 
we are concerned only with equatorial orbits. Conse¬ 
quently, we will set z max = 0. Instead of the parameters 
(^min; T’max) we use the semilatus rectum p and eccentric¬ 
ity e, defined by 


^min 
T max 


P 

1 + e ’ 



(18) 

(19) 


to indentify orbits. 


To compare orbits between different spacetimes we need 
a gauge invariant way of identifying orbits. A useful set 
of invariants is given by the orbital frequencies [22j , 


Q(j> 


2ir 


T r ’ 

<f> 

<£>’ 


( 20 ) 

( 21 ) 


where T r is Boyer-Lindquist coordinate time between two 
passes through periapsis, and the angular brackets (•) 
denote averaging with respect to proper time, 

(F) = lim J— f T F(t) dr. (22) 

T—>oo Z / J _j- 


For future reference we note that for periodic quantities 
this definition reduces to, 


(F) 


1 f Tr 

— / F(t) dr 

*r JO 

l f K 

- / K(A)EdA, 
tr Jo 


(23) 

(24) 


where T r is the proper time between two passes through 
periapsis, and A r is the corresponding amount of Mino 
time. 


The pair (fi^fi^) provides a gauge independent char¬ 
acterization of eccentric equatorial orbits. This charac¬ 
terization is almost (but not quite) unique. In a small 
region near the transition to plunging orbits, all orbits 
come in isofrequency pairs of distinct orbits with identi¬ 
cal frequencies [39]. 


D. Self-force corrected motion 


We now consider the situation that 0 < m « M(= 1). 
In that limit, we can study the effect of the small mass m 
on the gravitational field as a perturbation to the black 
hole geometry generated by the larger mass M, using the 
mass-ratio p = m/M as an expansion parameter. The 
corrections to the motion of the smaller object due to self¬ 
interactions can thus also be studied order-by-order in p. 
Detweiler and Whiting [U| showed that, to first-order 
in p, the mass m will follow a geodesic in the effective 
spacetime 

9iw = 3/ji/ + fyui/i (25) 

where g^ is the background Kerr geometry generated 
by the larger mass M, and h^ v is a certain smooth piece 
of the retarded metric perturbation generated by the 
smaller mass m. 

The regular piece of the metric perturbation, can be 
further split in a dissipative piece h^ lss , and a conserva¬ 
tive piece /i^ cons - By setting appropriate boundary con¬ 
ditions on the retarded metric perturbation, the effects 
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of the dissipative and conservative piece can be studied 
independently. The dissipative piece encodes all effects 
on the motion of m due to energy and angular momen¬ 
tum being carried away by gravitational waves causing 
the orbit to slowly decay and spiral towards the larger 
mass M. The conservative piece then encodes all other 
effects like changes in the orbital periods. 

Because the motion in the “conservative” effective space- 
time, 


~9™ ns = 9^ + Kr nS , (26) 

is free from dissipation, there exist bound eccentric equa¬ 
torial orbits in the effective spacetime with well-defined 
orbital frequencies 


Q(j> 


2n 



(27) 

(28) 


where the tildes denote quantities calculated in the effec¬ 
tive spacetime. 


E. Generalized redshift invariant 


To compare the conservative effect of the self-interaction 
between various computational approaches one needs a 
gauge invariant measure of this effect. In this paper we 
will calculate one such measure, the so-called (general¬ 
ized) redshift. This was first introduced for circular or¬ 
bits by Detweiler |I2!, and later generalized to eccentric 
orbits by Barack and Sago (221 [23]. Recent work by Le 
Tiec [H] has shown that the generalized redshift plays a 
key role in the ‘first law of mechanics’ for compact bina¬ 
ries, where it acts as a conjugate variable to the particle 
mass. In its generalized form the redshift invariant is 
defined as, 


U = 


& 


(29) 


periodic nature (or helical symmetry in the case of circu¬ 
lar orbits) of the spacetime [22] . Consequently, U viewed 
as a function of the orbital frequencies (fl r , f2^) provides 
a gauge invariant measure of the conservative motion. 
We can thus compare U in the background spacetime g^ 
and the conservative effective spacetime g™ ns (at fixed 
orbital frequency), 


a u(n ri n^) 


u(n r ,^)-u(n r M 

9 


(32) 


Following [23j we obtain U in terms of quantities that we 
can calculate directly for perturbation theory, 


U = 


9 \Tr T r J 

--^(fr-Tr) d r + OA) 


(33) 

(34) 

(35) 


From the normalization conditions for the 4-velocity 

n _ dx^ 

“ — dr > 


= -1, (36) 

one readily obtains 

^ = 1 -^ con W' + (lV), (37) 

and consequently 

< 38 » 

where h* u = = /i^ con Vu ! '. 


F. Radiation gauge 


In this paper, we will obtain the metric perturbation h^ 
in the so-called outgoing radiation gauge or ORG. This 
gauge is defined by the gauge conditions 


For circular orbits (and a helically symmetric choice of 
gauge) is independent of r and this definition reduces 
Detweiler’s original definition 

U=^- (30) 


For eccentric orbits we have, 





TV 

Tr 


(31) 


The redshift invariant U is invariant under gauge trans¬ 
formations that are asymptotically flat and preserve the 


= 0, and (39) 

K = ( 4 °) 

This gauge is naturally asymptotically flat, and the gauge 
conditions respect the periodic nature of orbits in the 
background spacetime. As such the ORG is a natural 
gauge to calculate AC/. 

A complicating factor of work in radiation gauges is that 
they are known [32] [33] to develop singularities in the 
presence of matter, and this singularity will generally ex¬ 
tend away from the worldline. This issue was studied in 
great detail by Pound et al. [33] . They found that radi¬ 
ation gauges fall in one of four classes. In the first two 
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classes, the gauge irregularity forms a string-like singu¬ 
larity extending along one of the principle null-directions, 
either towards infinity or the horizon. These are known 
as half-string gauges. In the so-called full-string gauge 
the singularity extends in both directions. Finally, one 
can stitch together two regular halves from the half-string 
gauges to form a gauge that is regular both at the hori¬ 
zon and infinity. The price paid is that this “no-string” 
gauge is discontinuous on a sphere containing the particle 
worldline. 

In this paper, we obtain the local metric perturbation 
generated by a particle by approaching the particle’s 
worldline either from the inside or outside in the half¬ 
string gauge that is regular on that side. The inside and 
outside limits are therefore obtained in different gauges. 
This will not cause problems as long as we remain aware 
that the limits of gauge dependent quantities may not 
agree. 


III. 

METHOD 

A. 

Strategy 


The core idea of our strategy is to separate the metric 
perturbation in modes depending on just one coordinate 
variable, and find these functions mode-by-mode. Such 
‘frequency domain’ methods have previously been used to 
great success to solve the linearized Einstein equation in 
a Lorenz gauge for small objects orbiting a Schwarzschild 
black hole [35] 3H 33 ■ Although this method has previ¬ 
ously been applied to scalar waves on a Kerr background 
[3DH35] . it cannot be applied directly to find perturba¬ 
tions to the metric on a Kerr background, because the 
linearized Einstein equation in Kerr is not separable. In¬ 
stead we have to work around this issue. 

Teukolsky famously showed that the linear equations for 
scalar fields with arbitrary spin-weight on a Kerr back¬ 
ground are separable (35). In particular, this holds for 
the equations for the Weyl scalars ipo and f/q, which have 
spin-weight +2 and —2, respectively. In a 1973 theorem 
[50] . Wald showed that the perturbations to the Weyl 
scalars ipo and ip 4 individually contain almost all infor¬ 
mation about the metric perturbations (up to global per¬ 
turbations of the mass and angular momentum). Cohen, 
Chrzanowski, and Kegeles (CCK) [MIP?7| developed a 
procedure to reconstruct vacuum metric perturbations 
in a radiation gauge from tpo or i/q. 

Our strategy will be first to solve the Teukolsky equation 
numerically by separation of variables, and then obtain 
the metric perturbations using the CCK reconstruction 
procedure. This technique was pioneered by the group 
of Friedman et al in Milwaukee [28113(71 [STl , applying it 
successfully to obtain the metric perturbations produced 
by a mass moving on a circular equatorial orbit around a 



Figure 1. Graphical representation of the method of extended 
homogeneous solutions. The field equations are solved in the 
interior and exterior vacuum regions in the frequency domain. 
To obtain the field at the particle the time domain solutions 
are constructed and analytically extended towards the world¬ 
line (either from the inside or outside). 

Kerr black hole [3D]. Here we are extending these results 
to eccentric equatorial orbits in Kerr. 

The extension to eccentric orbits brings new challenges. 
For one, the Teukolsky equation for the radial functions 
now has a non-zero source in a region extending from 7 ' m i n 
to r max . This causes a problem because the CCK recon¬ 
struction procedure is ill-defined for non-vacuum regions 
35. . Consequently, it cannot be applied mode-by-mode 
in the source region. 

We avoid this issue by applying the so-called “method 
of extended homogeneous solutions”. This method was 
originally introduced 135] to avoid poor convergence of 
the mode-sum at the particle for eccentric orbits due to 
the Gibbs phenomenon. The basic idea is as follows. A 
particle traveling on an eccentric orbit naturally splits 
the rf-plane into an interior and exterior region (see Fig. 
|T|). Outside of the region between r m j n and 7' max , the 
frequency domain equations are free of source terms and 
the field perturbation is given as a sum of homogeneous 
frequency modes with appropriate (retarded) boundary 
conditions at infinity (in the exterior) or at the horizon 
(in the interior). In the (1+1) time domain, the solution 
in the entire exterior/interior region is vacuum. So it can 
be obtained by analytically extending the frequency do¬ 
main solution in the appropriate vacuum region. This an¬ 
alytic extension can in fact be done on a mode-by-mode 
basis, yielding mode-by-mode contributions to the field 
perturbation at the particle. The sum of these contribu¬ 
tions (over the frequency modes) converges exponentially 

m- 

Schematically, our approach will be as follows. For each 
frequency mode we solve the Teukolsky equation to find 
the interior and exterior vacuum solutions for if) 4 . Al¬ 
though these will only represent the physical perturba¬ 
tion of 7/14 in the interior and exterior regions respec- 
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tively, these solutions will be valid vacuum solutions ev¬ 
erywhere. We then apply the CCK reconstruction pro¬ 
cedure to obtain the metric perturbation corresponding 
to these vacuum solutions. The time domain form of the 
frequency domain modes of the metric perturbation are 
then evaluated at each point along the orbit. The full 
metric perturbation is the sum of all these contributions. 


The plan for the rest of this section is as follows. In 
section IIIB[ we review the Teukolsky formalism for ob¬ 
taining i/> 4 . We pay particular attention to obtaining the 
asymptotic behaviour at infinity and the black hole hori¬ 
zon of the homogeneous solutions of the radial Teukol¬ 
sky equation. The first step in the CCK procedure is 
to obtain the so-called Hertz potential, which itself is 
a solution of the Teukolsky equation with spin s = 2. 
Section III C| leverages the asymptotic behaviour of the 
solutions to the Teukolsky equation to solve the equa¬ 
tions for the Hertz potential algebraically. The rest of 
the steps of the CCK procedure are discussed in section 
III Dt where we explicitly obtain an operator that pro¬ 
duces the (spherical harmonic) modes of h uu from the 
(spheroidal spin-weighted harmonic) modes of the Hertz 
potential. Section |IIIE| discusses the procedure for ob¬ 
taining the remaining pieces of the metric perturbation 
due to perturbations to the total mass and angular mo¬ 
mentum of the spacetime. This gives us the modes of the 
retarded metric perturbations, whose sum diverges at the 
particle. Section III F discusses the mode-sum regulariza¬ 
tion procedure for obtaining the regular “i?” part of the 
perturbations. 


B. Teukolsky Equation 


We here review the Teukolksy formalism focussing on the 
asymptotic properties of the solutions of the Teukolsky 
equation. For a more comprehensive review see [52l and 
the reference therein. Discussion of the static modes sim¬ 
ilar to ours can be found in [53] . 

A classical result from Teukolsky [19] IM] is that the equa¬ 
tions for fields <!> s with spin-weight s on a Kerr back¬ 
ground can be solved by separation of variables. Writ¬ 
ing 


= 


V / 27t, z 4 

L,m : n 


V s^lmn sRlmn( r ) sS t 


Imn 


(z)t 


lirrKp—ujjj 


(41) 

one finds that the radial mode functions sRi mn ( r ) and 
spheroidal mode functions 8 Si mn (z) satisfy individual 


2 Here we anticipate that we will deal with solutions with a discrete 

spectrum indexed by the integers m and n. In the more general 
case of fields with a continuous spectrum the sum over n would 
be replaced by an integral over lj. 


(uncoupled) ordinary differential equations, 
K 2 mn - 2 is(r - 1 )K„ 


{ A— 3 — | 

r A s+1 —) 

l dr' 

V dr/ 


A (42) 


+ 4 isLOmnT - 


s Arran ) s^Imn( r ) — s^ImnW) 


r d /. ,, d\ (m + sz . , 2 

— s(s — 1) + s A; mn j - sSi mn (z) = 0, 


(43) 


where 


Kmn = (r 2 + a 2 )u mn - am, (44) 

sArran = sA/ m n 4” CL LO rnn 2ui(XLO rnn . (45) 

The polar mode functions s S lmn {z ) are known as spin- 
weighted spheroidal harmonics. In the limit aoj mn —i 0, 
the eigenvalues reduce to 


sA/rran — sAmra — 4” 1) ^(s 4“ 1), 


(46) 


and equation (43) reduces to the familiar equation for 
spin-weighted spherical harmonics. In this paper we 
choose to normalize the spin-weighted spheroidal har¬ 
monics such that 


1 

sS lmui (z ) 2 = 1. (47) 

-i 


For our purpose, the most relevant examples of spin- 
weighted fields are the linear perturbations of the Weyl 
scalars 4> 2 = 4>o and d>_ 2 = p _4 t/q- The rest of this sec¬ 
tion will deal with the properties of the solutions of the 
radial Teukolsky equation (42) with |s| = 2. For reasons 
that will become apparent in Sec. |III C[ we will focus 
on determining the asymptotic behaviour of the homo¬ 
geneous solutions at infinity and near the horizon. 


We will distinguish three separate cases depending on 
whether uj mn and/or ma vanish. 


1. Homogeneous solutions (generic case uimn A 0 an d 
ma ^ 0) 


In the most generic case, where neither w m „ nor ma van¬ 
ish, the homogeneous solutions of the radial Teukolsky 
equation (42) are oscillatory. One straightforwardly es¬ 
tablishes that at infinity and near the horizon the homo¬ 
geneous solutions behave as [55] , 


sRimn{r)= sOtimn ^ 1 exp(-iu mn r*) r ->■ oo, (48) 
4 -sPimn r ~ {2s+1) exp (iuj mn r*) 
sRimn{r)= sAi mn A~ s exp(-ik mn r*) r r+, (49) 
+sB lmn exp(ik mn r*) 






















where k mn = u mn - and r* = r + r ^ r J r _ log ^ 


' 9 ' + is the tortoise coordinate. 


2r — log 

The positive and negative frequency modes (±w mn , or 
±k mn ) are naturally identified as representing incoming 
and outgoing waves at infinity and the horizon. The 
physical solution at infinity, sRtmni the solution with 
no waves coming in from infinity, implying s aiL n = 0. 
We further normalize this solution by setting s /3 + mn = 1. 
Similarly, the physical solution at the horizon, aRimn> 
is the solution with no waves coming “up” the black 
hole horizon, leading to the conditions = 0 and 

s Ajmn = 1- For later reference the complementary “un- 
physical” solutions at infinity and the horizon are de¬ 
noted a R* mn (defined by s /3f mn = 0 and s af mn = 1) and 
sRfmn (defined by „Af mn = 0 and s Bf mn = 1). 


Since (42) is a second order linear differential equation, 


its space of homogeneous solutions is two dimensional. 

I AL 

Consequently, any of the four solutions s Ri mn an d sRi mn 
can be written as a linear combination of two other solu¬ 
tions. In particular, we can write the “unphysical” solu¬ 
tions as a linear combination of the “physical” solutions 


Rf — 

L lmn 

s^lmn 

— R + 

Si^lmn s L lmn 

(50) 


> 



s ^lmn 


r /L — 

L lmn 

p+ 

sJri lmn 

— A + R~ 

S ' n ~lmn s x lmn 

(51) 


S B+ 


2. Axisymmetric static modes (o mn = 0 and ma = 0 ) 


The asymptotic behaviour of the homogeneous solutions 
changes in the special (static) case that oj mn vanishes. In 
that case (42) becomes 


0 = R"(x)+ {1 + °] {1+ , 2x) R'(x)+ 
w x(l+x) 

(a 2 — 2is{\ + 2x)a (l — s)(l +1 + s) 


4(x(l + x))‘ 


j(1 + x) 


)R(x), 


(52) 


These solutions have the asymptotic form, 
sRlmnfa) 

X s A i mn ((l — S )s + (2 — S) s -il(l + l)x 

+ (3-s) s _ 2 (/-1) 4 x 2 + 0(x 3 )) (54) 

+ s Bimn (divergent terms) 
as x —> 0, 


and 


sR‘lmni . X ) — 


„l — S 


(*Amn 


(20! 




sBlmn 


as X —> OO, 


mi-s) 

y/ns(l + f) s 1 


T+^)) + 


1 (x s A lmn 


(55) 


2 2/+l 


+ <^)) 


where X is an unspecified function of s and l that we will 
not need. Note the use of the Pochhammer symbol 

(x) n = r ^, + = x{x + 1) • • • (x + n - 1). (56) 

T(x) 


The physical solutions are identified by imposing regular¬ 
ity at the horizon and infinity. At first glance this may 
seem impossible, since both solutions in (53) appear di¬ 


vergent at the horizon. This is due to irregularity of the 
Kinnersley tetrad at the horizon, and the physical regu¬ 
larity condition is that A s R(r) should be smooth at the 
horizon [53]. We thus identify the physical solution at 
infinity, s R+ mn , as set by the conditions 


s A Ln = ° and 


(57) 

sBLn = (-l) S+1 2s(2«r s -^T, (58) 


(l + s)V 

and the physical solution at the horizon, s Rif mn , 
by 


as set 


s A i mn — (2k) 


and sB lmn = 0. (59) 


where x = r =Y ± , cr = - mo . and k = 

This equation has explicit analytic solutions. When the 
product ma also vanishes (i.e. for general static modes 
in Scliwarzschild (a = 0) or axisymmetric (to = 0) static 
modes in Kerr), eq. (|52| is solved by 


(53) 


sRlmn —s A lmn (x(l + X)) 2 Pf (1 + 2x) 

+ sB lmn (x(l + x)) 2 Qi (l + 2x), 


where P™(x) and Q™(x) are associated Legendre P and 
Q functions. 


3. Static modes (uimn = 0 and ma ^ 0) 


We finally consider the very special case that w mn = 0 
while both a and to are non-zero. For modes sourced 
by a particle moving on an equatorial orbit this can only 
happen if an integer combination of the radial and az¬ 
imuthal frequencies Q r and f2^ vanish. Such special or¬ 
bits are known as r</>-resonances, which can cause a co¬ 
herent build up of linear momentum flux to infinity |56j , 
but should not have a direct impact on the local dynam¬ 
ics of the binary system at leading order in /. i. 
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The most general solution in the case that ma ^ 0 is 4- Inhomogeneous modes 

given by, 


sRimn s A , 


lm 


l + x 


(x(l + x))~ 


x — S' 1 + l — s; 1 + s — ia\ —x) 


sB irn 


l+x 


x 2 Fi ( l + s; 1 + l + s; 1 — s + ia; — x ), 


(60) 


where 2 -P 1 a hypergeometric function. 

Near infinity and the horizon this solution has the fol¬ 
lowing asymptotic form 

s R lmn( X ) = Z~ S_lf (sAm + °( X )) 

+ x 7,2 ( s B lm + O(x)) (61) 

as x —» 0, 


and 


The solution of the inhomogeneous Teukolsky equation, 
sourced by a particle moving on an equatorial orbit, in 
the vacuum regions inside (—) and outside (+) the orbit 
can be written as a linear combination of the homoge¬ 
neous solutions, 

^=iE Z Ln-2 R Ln(r) 

v l,m,n 

(69) 

The coefficients Z^ u can be determined using variations 
of parameters, 


-± = r max - 2 ^ n (r)- 2 r tmn (r) 

W L* WUR+ mn ,- 2 RT m J(r) ’ 


(70) 


where W[fi, f 2 \ denotes the Wronskian of the solutions 
/1 and f- 2 . An explicit expression for the “source” term 
- 2 TlmrSr) is s iven in E 3 - 


s^lmn 


(x) = x l -°( 

+ sB { 


(1 + l + s)l—s 

(1 - S ~ l<T)i+ 
(1 + 1 — s)z+£ 


(1 + s -(- ia)i- 
+ x 1 s 1 (X s A lm + 
as x —> 00 , 


Oil) 


Y s B lm 


(62) 


O(D) 


where X and Y are unspecified functions of s, l, and a 
that we will not need. 

Again imposing regularity, we identify the physical solu¬ 
tions at infinity s R+ mn by 


4+ 

s ^lmn 


(0 A _ s (l~s)\(l-S-i<T)l +s 
(l + s)!(l + s + icr)i- s 


s Bf 

* Lmn 


= 1, 


and 


(63) 

(64) 


and the physical solution at the horizon, s R i mn , by 


s A imn = ( 2k ) S and sB lmn = 0. (65) 

Moreover, we will need the irregular solution at the hori¬ 
zon s Rf mn defined by 

sAf mn = 0 and a Bf mn = ( 2 k)~ s , (66) 


which can be written as a linear combination of the phys¬ 
ical solutions at the horizon and infinity, 


TD~^~ _ A + 

r>-f _ s± 'lmn 

sJr ^lmn 0 + 


TD— 

mn *' lmn 


,b 


sR lmn A 


Imn 

(l - s)!(l - s - ia)i +s 
( l + s)!(l + s + i<y)i~ s 


Ji, 


(67) 

( 68 ) 


C. Hertz Potential 


The reconstruction of the metric perturbations in the 
CCK formalism is defined in terms of the so-called Hertz 
potential. In vacuum regions, the Hertz potential for 
metric perturbations in the outgoing radiation gauge 
+ORG satisfies two conditions [55]. First, 'I’org satis¬ 
fies the Teukolksy equation for s = 2 fields in vacuum. 
Second, +org satisfies a fourth order differential equa¬ 
tion with tpi acting as a source term, 

^A 2 (Dj) 4 A 2 4- 0 RG = p-V4, (71) 

where D\ = d r — b + a )9t+ad^ , anc j t ]. le over j jar denotes 
complex conjugation. 

The first condition implies that +org in the interior 
and exterior vacuum regions can be decomposed in spin- 
weighted spheroidal harmonics, 




± _ 
ORG ~ 



X! ^tnnzRtrnnir) 2S lmn {z)e i( ‘ 

l,m,n 


... (72) 

Furthermore, the linear operator on the left hand side in 
( [7l] ) neatly separates over the spheroidal modes, with the 
operator acting on the radial modes given by 


Bmn. — d r + 




i(r 2 + a 2 ) - ma _ n 

— —m—n • 


A 


(73) 


Consequently, <zu can be solved mode-by-mode. In par¬ 
ticular, we observe that the resulting operator acting on 
the radial modes is a linear operator that maps solutions 
of the radial (vacuum) Teukolsky equation with s = 2 
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into solutions of the radial (vacuum) Teukolsky equation 
with s = —2. Consequently, if we fix a basis on the (2- 
dimensional) space of solutions of the (vacuum) Teukol¬ 
sky equation (as we did in section III BI, the action of 
this operator is given by a 2-by-2 matrix, which can eas¬ 
ily be inverted. The basis in section [III B| was chosen to 
ensure that this matrix is diagonal, simplifying the in¬ 
version. Moreover, the action of the radial operator can 
be determined at any radius r. In particular, we can 
determine the action (and its inverse) in the asymptotic 
regions at infinity and near the horizon, requiring only 
the asymptotic behaviour of the homogeneous solutions 
as obtained in the previous section. 

Ori [23] performed this analysis for the relation between 
(generic) modes of the incoming radiation gauge (IRG) 
Hertz potential and i/jq. We here extend that analysis to 
the ORG Hertz potenital, 'S’org, and ip 4 , also including 
all the static modes. 

Expanding the left hand side of © using (|72|) we find 


— A 2 (Z?J) 4 A 2 'I' 


32 


ORG 


E 

l.m.n 




A 2 (Z) t ) 4 A 2 

\ —m—n J 


X 2 RLn( r ) 2<S) 


Imn 


32 

( Z ) 


(74) 


Relabelling (to, n) to (—to, — n) and using the identi¬ 
ties sS lmn {z) = (—1 ) s+m _ s S)_ m _ n (U) and 
sRi- m - n , this becomes 


Teukolsky equation with the opposite spin —s, as can eas¬ 
ily be verified by inserting A s f(r) in (421 with spin — s. 


For generic (uj ^ 0) modes, the asymptotic behaviour 
tells us that this relation exchanges physical (retarded) 
for unphysical (advanced) boundary conditions, 

A S sRLn= sRfmn- (79) 

Using this relation Eqs. (76) and ([77f become 


-2 Rhnn( r ) — 

E§± A 2 (Dt ) 4 nRf (r) 

22 l—m—n K-^mnJ ~^ L l—m—n\ ) 


(80) 


and 


Zlmn(Dmn) A - 2 Rl mn { r ) ~ 

(-!)"*,T,± 


32 


-2 Ri_ m _ n (r). 


(81) 


If we insert the asymptotic expansion (48) into (80) we 
find 


Z Ln r3 exp(fw m „r*) = 




32 


^(U m _ n 16 u;V exp(iw m „r*), 


(82) 


= for the solution at infinity. This can easily be solved for 


\I/+ 

Imn ’ 


^A 2 (^) 4 A 2 T± i?G 


E< 

l,m,n 


-iy 


, ^tm- u . 

32 


*Ln = (“I) 


y+ 

£+mn ^Imn 




(83) 


A 2 (£)t l J 4 A 2 2 Rf m 


,(r) ^ 2 S lmn {z)e i( - rn ' 1, (75) where we used that Zi mn = (-l) l Zi_ m _ n . 


Comparing with the expansion (691 of the right hand side 


of (71), and noting the orthogonality of the spheroidal 


and Fourier modes, we obtain, 


J lmn ' 


1 Iran M = 

(-D 


32 


^_ m _ n A 2 (0 4 A 2 2 i?L(r). 


L lmn V 

(~i) r 

32 


^ l-m-nP lmn ‘ 2 Rlmn( r )’ 


(76) 


By the Teukolsky-Starobinsky identities (see e.g. 
section 81) this is equivalent to 

Z Ln(D mn ) A -2RLn( r ) = 


(77) 


However, for the physical mode at the horizon if one in¬ 
serts the asymptotic expansion near the horizon ( |49| ) into 
(80), one finds that the leading terms cancels, yielding no 


useful information. Using (81) instead one finds, 

^mn( 2K ) 4 (*( a + 2u -»») ~ 2) 4 exp (~ik mn r*) = 

(-l) m -_ ( 84 ) 

—^—^i- m - n Pimn exp (-ik mnT ) 1 


which we can solve for T 


Zmn’ 


907 “ 

Kmn = (-l) z+m —^(2 K ) 4 (*(a + 2u mn ) 2 ) 4 . (85) 

Plmuj 


where 

Plmn ~ (( —2 ^Imn T 2) T 4m.ttCU mn 4d CU mn ) 
X ( —2 ^Imn T 36TOUCU mrl 36u 
T (2 —2^imn T 3)(9 Qa“uj rnn 48TOOW mr 
+ 144^(1 - a 2 ). 


(78) 


We now observe that if f[r) solves the radial Teukolsky 
equation (42) with spin s, then A s f(r) solves the radial 


For the static case to = 0, but ma ^ 0, one easily verifies 
that 

a s p+ _ {l ~ s)\(l - s ~ ia)i +s n+ 


(/ + s)!(l + s + icr)i— 8 


-sRtmu, (86) 


and 


A* sRlmn — -sRlmn- 


(87) 
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Repeating the analysis above we find 

32 


\l>+ = (-1)' 

lmn v / 


(icr — 2)4 


7-1- 
^Imn ’ 


*Zmn = f- 


/ i\Z+mqo ^) 4 7- 

( } lmn ' 


Similarly, for ui = ma = 0, 


^ s-^lmn ~ y _|_ s ), an( ^ 

AS D- _ + S ) ! D- 

^ sJx lmn ^ _ g ^j — sIx lmn > 


( 88 ) 

(89) 

(90) 

(91) 


leading to 






lmn 


= (-1) 


Z+m+lno 7+ 


32 


((i -1)4) 


7 - 

2 imn' 


(92) 

(93) 


This brings us to the main result from this section: a set 
of algebraic expressions for the asymptotic amplitudes 
T f mn of the outgoing radiation gauge Hertz potential in 
terms of the asymptotic amplitudes Z^ nn of ^> 4 . In sum¬ 
mary, these expression are, 


1 

( ( 1 \l+m 2 7 + 

\ L > uj^ nn Z lmn 

for 

Mmn 


o, 



^lmn ~ \ 

' ( 1 y+m 32 ry-\- 

I V L > (icr— 2)4 n lmn 

for 

Mmn 

= 

0 but 

ma 

7 ^ 0 . 

\ 

{(-iy+^32 Z+ mn 

for 

^mn 

= 

ma = 

0 , 



U-l)l+™J^( 2 K y(i{o + 2 u mn ) 

- 2)4 Zr mn for 

^mn 

7 ^ 

0 , 



^Tmu = 


for 

^mn 

= 

0 but 

ma 

7^0 

\ 

( 1 \l+m 32 7 - 

l 1 " ((Z-l ) 4 ) 2 mrrm 

for 

^mn 

= 

ma = 

0 . 



(94) 


(95) 


In addition, by combining (150)), (|51l), (67), (79), ( 86 ), (87), (90), and (91), we can obtain the s = 2 homogeneous 


solutions of the radial Teukolsky equation needed in (72) explicitly in terms of the s = — 2 homogeneous solutions, 


2 R LW 


2 


R lmn( r ) 


1 

A 2 


V -2 <*i 


"=~ - 2 RTmnir) - -2RLn(r) 


(1 —2 + 8<T)i + 2 


A 2 (Z-l) 4 (l+2-i<T ),_ 2 


-2 RLu( r ) 




1 

A 2 


-2 RLn( r ) - 


-1 R lr 


A 2 1 - 
0 - 1 ) 

A 2 


u . 

( D+ I (^~ 1 )4(l+2+io-)|-2 D— \ 

y ~‘^ rt lmn ' (1—2—icr)z_|_2 ~‘^ Ix lmn J 


»(0 


~2‘Rlr. 


for uj rnn 7 ^ 0 , 

for w mn = 0 but ma 7 ^ 0 . 

for oa mn = ma = 0 , 

for uz mn 7 ^ 0 , 

for w mn = 0 but ma 7 ^ 0 , 

for w mn = ma = 0 . 


(96) 


(97) 


Consequently, once we have obtained the mode expansion 
(69) for - 04 , we can use these algebraic relations to ob¬ 
tain 'i’oRG (in the interior and exterior vacuum regions) 
without solving any further differential equations. 


The CCK procedure obtains the remaining non-zero com¬ 
ponents of the metric perturbation by acting on 4 >org 
with certain second order linear differential operators (see 
e.g. [58] V With our sign conventions these are given by 


D. Metric Reconstruction 

The next step in the CCK procedure is to produce the re¬ 
constructed metric components from the Hertz potential. 
These will be produced in an outgoing radiation gauge 
defined by 

(98) 

(99) 



hu = 

e l e l h/iv = 'Hll'&ORG + C.C., 

(100) 


h 13 = 

e i e 3 ^ = 'Hu'&ORG, and 

(101) 


h-33 = 

e 3 e 3 = ii33^0RG, 

(102) 

with 




H 11 

II 

1 

^1 

1 

3 a — j3 + 57t) (5 — 4a + 7r) , 

(103) 

R-13 

= >-{(* 

— 3a + fl + 57t + f ) (A + /a — 

4 7 ) (104) 


(A + 5/Li — /j, — 3y — 7 ) (5 — 4a + 7 r) j, and 


/l2a — ^2 — 0? 

/I34 = 6364/2^1/ — 0. 
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%33 = P 


-4 


(A + 5^ - 3y + 7 ) + fj, — 47 ^ . (105) 


Here 5 = e^d^, A = and the remaining Greek 

letters are the usual spin-coefficients in the Newman- 
Penrose formalism. (Explicit expressions consistent with 
our conventions are given in El)- 

Our goal in the rest of this section is to construct h uu 
from <F org, 

h± u = lim u a u b h ab = lim Huu^org + c - c - ( 106 ) 


with 


H uu = m 1 w 1 'Hh + 2it i M J ‘Hi 3 + u 6 u' i 'H^, 


1„,3- 


,,3„,3<] 


(107) 


and where the limit x —► symbolizes approaching the 

particle worldline from either the exterior (+) or interior 
( —) vacuum region. As written here Eq. (106) makes lit¬ 
tle sense, since the particle four-velocity u M is not a field, 
but a quantity defined only on the particle worldline. To 
make sense of the limiting procedure we must extend 
to a suitably smooth held defined in a neighbourhood of 
the worldline. Following the literature ei mi, we choose 
a “rigid” extension, i.e. we choose to extend the four- 
velocity iA such that 


u^ih, r,z,<t>) = u>*{to,rQ,zo,(l>o), 


(108) 


using that the particle crosses each t 0 -slice at a unique 
point (fo, zq, 


A second issue with (106 1 is that the limit towards the 


particle worldline will obviously diverge. However, to 
obtain A U we only need the regular “R” part of h ab . 
To obtain the regular part, we will utilize the mode-sum 
formalism [6D - I63] . This requires that we expand h uu in 
spherical harmonics, whereas we obtain ^org expanded 
in spin-weighted spheroidal harmonics. The remainder of 
this section will be devoted to explicitly evaluating the 
action of 'H uu and “re-expanding” the result in spherical 
harmonics to obtain the “Z-modes”. 

We start by inserting the expansion of the Hertz potential 
([721) in (fl06|), 


h iu = iim , - 7 == 

x—>Xq V Z7T 


E 


R Ln( T ) 


l,m,n 


x 2 S 1 „(z)e I W-“”” f )}+c.c. 


(109) 


The first step towards re-expanding the result in spheri¬ 
cal harmonics is to expand the spin-weighted spheroidal 
harmonics s S lmn (z) in spin-weighted spherical harmon¬ 
ics s Y lm (z), 

sSi imn (z) = E ](sbmn)il sYl 2 m( Z )- (HO) 


Hughes [64] has described a robust way of numerically 
obtaining the transformation matrix (, sb mn )i 1 - 


After expanding to spherical harmonics, we can evaluate 
the action of the linear operator H. U u- Observing that 
derivatives with respect to z can be rewritten in terms 
of the spin-weight lowering operator = y /1 — z 2 {d z + 
1372 ^ + jz^ 2 ), the result has the schematic form 

ht= lim^dW—"^ E C hmnsi(r,z) 

x ^ x o m,n h ,l 2 ,i,s (HI) 

2Rt 1 ’^ n (r){2bmn) l l 1 2 S Y h m( z ) + C-C., 

where C is a numerical coefficient depending on the field 
point (r, z), the particle position (ro, Zq ), and the indices 
( li,m , n, s, i). 

Next, we expand the spin-weighted spherical harmonics 
s Y lm (z ) to ordinary spherical harmonics by using the fol¬ 
lowing identities 


2m.Ah y ( \ 

tfU(*) = E-T 2 


1 Y hm {z) = E 


— z- 

Im/il 


W 2 Y hm (z) 


t y/ih-Wi + yvr^’ 


0 Y lim (z) = E 


Y hm (z) 


i 2 \J(h — l)h(h + l)(h + 2 ) 

where (using Wigner 3j-notation) 

2 m A\i = (-ir\H^h + mi 2 + i) 

2/i I 2 \ / 2 ^1 ^ 


( 112 ) 

(113) 

(114) 


y 0 m — mJ \2 —2 0 J ’ 

(115) 

lm A\\ =(-l) m+ V2(Zi-l)(Zi+2)(2Z 1 + l) 


x \J 2 i 2 + 1 ^ 


1 h h \ h h 

0 m —m I \ 1 —1 0 ) ’ 


and 


0 m A\l = y '(l 1 -l)l 1 (l 1 + l)(h + 2)6 hl2 . 


(116) 


(117) 


The normalization of the coefficients sm A 1 has been cho- 
sen as to cancel the Z 1 dependence of the coefficients 
Ci im nsi■ The result is 

htu = lim E E C mnai (r,z) 

x ~>- x 0 rn,n i,s 

X E Kmn*K™n{r){2b mn )\l (H8) 

hihih 

x sm A\lY l 3 m {z) + c.c ., 

where C is now a different coefficient that (among other 
things) still depends on the field point z. Observing the 
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identities 


C mnsi 5 z) 

— ( 1) ^'—m—nsiiX') ^)i 

(119) 

limn 

= (-1 ) l K- m - n , 

( 120 ) 


= 2R^ ) m _ n (r), 

( 121 ) 

(2 b m n)i 2 

= (-l) ll+h (2 b-m- n )\l. 

( 122 ) 


sm A\l = (-1 )«+fa+*3»-m 4 ja ) ( 123 ) 

Y hm (z) = (-l) m Y la _ m (z), (124) 


and relabelling appropriately, we can evaluate the com¬ 
plex conjugate terms, 


htu = lim £ £ (C mn 8 i (r, *) + (-l)^ +m C mns! ,(r, -*))*± ron (r)( 2 6 ro „)^ m ^ h 3 ™W ( 125 ) 


0 m,n 


h,h,l3 


= lim .J2 eiim *~“ mnt) E 3+ ”* F i3m (z), 


(126) 


0 m,n 


i,s 

h,h,l3 


r 


where in the last line we have redefined the coefficients 
C, which are now only a function of z 2 . 

To eliminate the dependence on the field point z in the 
coefficients we utilize the freedom granted by the pres¬ 
ence of the limit towards the worldline. Since we are 
only interested in the limiting value at the worldline, we 
are in principle free to multiply each term in the sum in 
(1261 by a function /(r, z) that smoothly approac hes 1 
at the worldlinej^] For each term in the expansion (126) 
we choose / such that 


utilizing the identity 


^ Y hm (z) = ]T Y l2mn (z), (128) 


with 


lm B h = (_ 1 )m+b+ 1 (^ _; 2 ) 


and 


l\ + I 2 + 1 f 1 l\ l 

0 m —m 


(129) 


\ Cmnsi(j’ Oj^o) 

ffaz) = ^- 7 - 

Cmnsi y ■> Z ) 


(127) 


That is, we keep only the leading constant term of each 
of the coefficients C. 

We finally rid ourselves of the final dependence on z by 


imgii _ /lmgii'jJ 

fc 2 V fc 2 / 


i-2 




k =1 


(130) 


57 ^(n 1 ^) 1 "^- 1 - ( isi ) 


This allows us to completely evaluate the limit and obtain 
the /-modes of h uu , 


h± = _ 


E^^^olE E Ci 3ronsi (r-o^o 2 )^^„2<^(r 0 )( 2 6 mri )^”i4; 


l-(-l) l 3 + TT 


m,n 

1 


/1 ,Z 2 , 

^3 ?^4 


= E 7 =Ee l( ^°-“™ to) Ec si (m,n,ro) E 2 iS(r 0 )( 2 &™){”"Mi a *W(0) 

l * ^ m,n i,s h,h 

= V /i*’* 


t3 l ilY Um (z 0 ) (132) 

(133) 

(134) 


3 Although any choice is allowed a priori, the price we pay is that 
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where in the second line we utilized that kj 4 , m (zo) van¬ 
ishes on the equator for odd values of U + m, and conse¬ 
quently that 


E 


1 —( —l) i 3+ T1 


n i 

E 3 Y Um(V) = y h m(0). (135) 


Equation (133) thus gives the explicit expression for the 
i-modes of h uu that we were after. Following the above 
steps the coefficients C S i(m,n,r 0 ) can be computed ex¬ 
plicitly, with their final form given in Table [I] The proce¬ 
dure in this section was purposely written more general 
than strictly necessary for the calculation of h uu . In prin¬ 
ciple, the procedure described here can be applied to any 
quantity constructed from the metric and its derivatives, 
such as the self-force. This will just yield a different set 
of coefficients C S i(m,n,ro)- For more general quantities, 
the expansion coefficients will not disappear from 

the final expression. 


E. Completion 


Although the Weyl scalars ipo and ip 4 contain most of 
the information about linear metric perturbations, they 
cannot contain all information. The simplest coun¬ 
terexamples are perturbations of the background Kerr 
metric within the Kerr family. These necessarily have 
ipo = ip 4 = 0, and can thus not be distinguished from the 
zero perturbation on the basis of ip o and ip 4 . The CCK 
procedure can thus at best return a representative of the 
equivalence class of metric perturbations with the same 
ipo and ip 4 . In general, we can write 

h)™ = h c ™ + p , (136) 

where h is the metric perturbation obtained by con¬ 
structing ip^ from and then applying the CCK pro¬ 
cedure. The procedure described in the previous sections 
will produce (the regular part of) h^ K produced by a 
particle moving on an equatorial orbit. However, to cal¬ 
culate A U, we need h^ 1 ; i.e. we need to find /i“° mp . This 
is known as the “completion problem”. 

The completion problem is made tractable by a very use¬ 
ful theorem due to Wald m- Wald showed that, up to 
gauge terms, the completion part in any vacuum region 
is given by, 

^° mp = + cjh ^ + cch + c NUT h^y T , (137) 


we have to account for the choices made in the calculation of 
the singular part in Sec. 0 In practice, we want to make 
our choice such that the leading order divergent structure of the 
singular field encoded in the regularization parameters is unaf¬ 
fected. 


where the metric perturbations are obtained by 
embedding the background Kerr metric in the four- 
dimensional family of vacuum Plebanski-Demiariski met¬ 
rics [651 and varying with respect to the mass M, angular 
momentum J = Ma , C-metric acceleration ac, or NUT 
charge Qnut, and the c» are constants. 


This reduces the completion problem to finding four num¬ 
bers in each of the two vacuum regions in our problem. 
The perturbations /lA, and h^y T have conical singular¬ 
ities extending from the horizon to infinity. Since the 
metric perturbation in the interior vacuum region must 
be regular at the horizon and the perturbation in the ex¬ 
terior vacuum region must be regular at infinity, we must 
have Cq = c^[, r = 0. 

Requiring that the perturbed spacetime has the correct 
ADM mass and angular momentum fixes the parameters 
in the exterior vacuum region, c\j = m£ and cj" = mC. 
We obtain 


, /£ r, r 2 ((r + a 2 )£ — aC) „ 

+ = 2m( 7 ^ U J - L dr 2 

+ q((r + 2)Z-Mr + l)g) d0 ,_2g dtd 

r r 


(138) 


By requiring continuity of certain gauge invariant quan¬ 
tities away from the equatorial plane, Merlin et al. j.3fil 
have shown that for particles moving on equatorial orbits 
c m = c j = 6 - 

Wald’s theorem only fixes the completion part up to 
gauge modes. Recall that to calculate U we must be 
in a gauge that is asymptotically flat and respects the 
periodic structure of the orbit. The two components of 
the metric perturbation in the exterior vacuum region, 
h™ K ’ + and naturally satisfy these conditions. 

On the other hand, the reconstructed metric perturba¬ 
tion in the interior vacuum region h^ K, ~ has (gauge) 
string like singularity extending towards infinity [34] . We 
must therefore construct A U from the values of h uu ob¬ 
tained from the exterior vacuum region. 


F. Regularization 


The procedure above calculates constructed from 
the retarded solution to the linearized Einstein equation. 
This quantity obviously diverges at the location of the 
particle. This manifests itself in (133) as the outer sum 
over l diverging, while the individual Trnodes, h l uu are all 
finite. 


A key result of the self-force program is that the retard 
metric perturbation can be split into a smooth regular 
“R” piece, and a singular “S” piece, in such a way that 
the contribution of the singular “S” piece to the equations 
of motion of the particle vanish (see [4j and the references 
therein). Moreover, this split can be arranged in such a 
way that the regular “R” piece is a vacuum solution to 
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i s 


C is (m,n,ro) 


~2 —I— T~ 

r 0 u u 


0 0 
0 1 
0 2 


1 1 
1 2 

2 2 


—2aro(r 0 CJ mn + i)u x u 1 — \p2 [ir 0 K mn + 2(r% — a 2 )) m 1 u' 
a 2 rouJmn(roLOmn + 2i)« 1 u 1 + y/2auj mn ( ir 0 Km n + 2{rl - a 2 )) u 1 u 3 

/ A _ j Tf 1 \ 

+ (Aq(-+ ir 0 u> mn - 2) - - (K mn - 2i(r 0 - 1)) (A' m „ - 4i(r 0 - 1)) J 


3 3 
U U 


ro 

— V^Aoro ^ 1 ^ 3 


/ A \ 

V / 2aA 0 r 0 a; mrl -u 1 u 3 +^ — + A 0 ( iK mn + 4(r 0 — 1 ))Jm 3 u 3 


±o 3 3 

— u u 


Table I. This table shows all the non-zero coefficients C s ;(m, n, ro) of h uu in Eq. (1331. Here ro is the radial coordinate at the 
particle worldline, and Ao is A evaluated at r = ro. The u l are the tetrad components of the four-velocity of the particle. 


the linearized Einstein equation and the particle follows 
a geodesic effective space time 


~ lR 

9[iv = 9 + i^tiv 

(139) 

Consequently, as mentioned in Sec. |IID[ A U is defined 
with respect to this regular effective spacetime. Hence 
we want to calculate 

h uu(xo{\)) = lim ■ (h^(x) - hl u (x)). 

x^xo{X) 

(140) 

Both and h^ u diverge at the particle, 

individual Z-mode contributions to and h% v 

at the particle one can write following [61| , 

Since the 
are finite 

huu( X o) = lim ( h uu( x ) - h tu( X )) 

X—>Xq 

(141) 

OO 

= E ( h uu’ l ( x o) - hH(x o)) . 
1=0 

(142) 

A key idea behind the mode-sum regularization scheme 
is that the Z-modes of the singular field can be written as 

«-* + I + 1 +o( (i + V 

(143) 

In general this expansion depends on the chosen gauge 
and extension of vA and more generally the extension of 
the individual terms of h uu . It is most easily calculated 
in the Lorenz gauge. Working in that gauge (and the 
same “rigid” u^ is constant extension as used in Sec. 
HID), one can show [BD] that for a particle moving on an 
equatorial orbit around a Kerr black hole, 

1/if £2+a2+ TT ) 
j E+rg+a 2 + — 

(B=) g 0 L ° r = , - ° 3 , 

7 r v / £ 2 +r 2 + a2 + ^ 

(144) 

(<7 =) = 0, 

(145) 

OO ft 

(D=) J2h S uu-£°~ tEt =°> 

1=0 1 + 2 

(146) 


where the letters in parentheses are the traditional names 
for these “regularization parameters” used in [SUED [62]. 
Vanishing of the C and D parameters means that the 
mode-sum for the regular part of h uu simply becomes (in 
Lorenz gauge), 

OO 

C = E i h uu ’ 1 - £o) ■ (147) 

1=0 

However, the metric reconstructed using the CCK proce¬ 
dure is obtained in an outgoing radiation gauge. A main 
obstruction in radiation gauge is that in general the full 
retarded metric perturbation is also singular away from 
the particle with a string like singularity radiating from 
the particle in one of the principle null directions. 

An in depth analysis of this problem was given by Pound 
et al. in [53] , They calculate the radiation gauge reg¬ 
ularization parameters of the self-force by finding a lo¬ 
cal gauge transformation that transforms the radiation 
gauge to a gauge that locally near the particle resembles 
the Lorenz gauge. This “locally Lorenz gauge” falls in the 
class of gauges related to Lorenz gauge by a continuous 
gauge transformation, for which one can show that one 
can use the Lorenz form of the mode-sum formula. Pound 
et al. find that for self-force in one of the “half-string” 
radiation gauges (regular in either the interior or exterior 
region) the A , B , and C parameters are unchanged, but 
the D parameter gains a finite correction. 

Using the results from | f3'4] it is elementary to extend their 
analysis to the regularization parameters for (h uu ) in ra¬ 
diation gauge. Let be the gauge vector that transforms 
the radiation gauge to a locally Lorenz gauge. Under its 
action the change in h uu is 



(148) 

= 2 

(149) 

= 2«"V,K^) - 2^^V m wU 

(150) 


The first term is a proper time derivative of a spacetime 
scalar. Consequently, it will vanish after averaging over 
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the orbit. In [33] it is shown that 

^ = O (logs), and (151) 

u> i V ll u v = <D(s), (152) 


From a, p, and e we calculate the values of the (specific) 
energy (£), (specific) angular momentum (£), orbital fre¬ 
quencies (O r and fi^), and periods (T r , T r , and A r ) using 
the analytical formulas from EM- 


where s is the separation from the worldline. Conse¬ 
quently, the second term in (150) vanishes as one takes 
the limit to the worldline. We thus find that (h uu ) is in 
fact invariant under the gauge transformation defined by 
Therefore, we find that the correction to the regular¬ 
ization parameters for (h uu ) in radiation gauge must be 
zero. 


The mode-sum formula (h uu ) in radiation gauge is thus, 


OO 

(^ Rad!± ) - E «C t>RacM,± > - <£o° r » ■ 
1=0 


(153) 


We can apply this formula to numerically obtain the reg¬ 
ular field from our reconstructed metric. 


IV. NUMERICAL IMPLEMENTATION 


We have developed a numerical code in Mathematicci 
which implements the above method for a mass follow¬ 
ing an equatorial geodesic around a Kerr black hole. The 
code uses the following steps: 


1. Calculate the geodesic sourcing the g ravita tional 
perturbation. (To be discussed in Sec. IV A| ) 

2. Solve the Teukolsky equation to obtain the linear 
perturbation to ij} 4 . (To be discussed in Sec. IV B ) 

3. Reconstruct the linear (retarded) metric perturba¬ 
tion and h~ u . 

4. Regularize /i„„ using the mode-sum formula. (See 

Sec. |5 If1 ) 

5. Estimate the contribution from the high l tail. (To 
be discussed in Sec. IV D ) 

6 . Add the contribution from the completion. (See 
Sec. |5lE|) 


Step 1 only needs to be done once, and is done in a 
matter of seconds. Steps 2 and 3 need to be iterated 
over all modes. Section ITV Cl discusses the criteria used 
to truncate the infinite sums over modes. Once this loop 
is complete the remaining steps can be done in one go. 


A. Geodesics 


We start our calculation by specifying a spin a for the 
central black hole, and an eccentricity e and semilatusrec- 
tum p, which define an equatorial bound geodesic orbit. 


A key observation used throughout the code is that all 
the quantities that we are interested in are smooth peri¬ 
odic functions along the orbit. As a consequence, numeri¬ 
cal calculations of orbital integrals and averages converge 
exponentially with the number of sampling points, if one 
uses an even sampling. To make maximal use of this, we 
fix an evenly spaced (in Mino time) grid of N points along 
the orbit, and precomputed as many quantities along the 
orbit as possible. The number of points N will determine 
the precision of the orbital integrals used to obtain the 
orbital averages and the coefficients Z^ nn of the inhomo¬ 
geneous solutions of the Teukolsky equation. In practice 
n will range from 1 for circular orbits to about 400 for the 
most eccentric orbits (e = 0.4) that we calculate here. 

We start by computing the geodesics. Fujita and Hikida 
m . 1 have derived analytic solutions for bound geodesics 
around a Kerr black hole in terms of elliptic functions. 
Evaluation of elliptic functions is relatively costly. To 
save computation time we evaluate the analytic solutions 
(and their first derivatives) on our orbital grid and save 
these as arrays for future use. 


Further savings can be made by precomputing some func¬ 
tions along the orbit that will be used frequently. For 
example, the coefficients in Table [I] can be split in terms 
independent of Z, m, and n, which are functions along the 
orbit that can be computed along the orbit. To avoid 
recalculating these functions for each mode, we evalu¬ 
ate them once on the grid at this stage of the compu¬ 
tation and cache their values. A similar procedure can 
be used for the coefficients appearing in the source term 
- 2 T lmn (r) in m (see |57|). 


B. MST method 


We solve the radial Teukolsky equation using the semi- 
analytical method developed by Mano, Suzuki, and Taka- 
sugi (MST) [671168] , see [12] for a review. In this method, 
solutions to the homogeneous radial Teukolsky equation 
are obtained as a series expansion in certain hypergeo¬ 
metric functions. Details of our numerical implementa¬ 
tion can be found in HU, which refines numerical meth¬ 
ods developed in [53HTT] . 

When presented with the orbital data calculated in the 
first step, the code detailed in m will return the homo¬ 
geneous solutions — 2 Rimn and their first radial deriva¬ 
tives evaluated at the grid points along the orbit. In 
addition the code returns values for the coefficients of 
the inhomogeneous solutions Zjj nn , and the coefficients 
of the asymptotic expansion at infinity and the horizon 
of the homogeneous solutions, -2 Af mn , - 2 -B,mn» -2 
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and —2/3,1 


We can then use Eqs. (941,(95),(96), and (97) to con¬ 
struct the value of the modes of the Hertz potential at 
each grid point along the orbit. The first radial deriva¬ 
tives are obtained from the radial derivatives of Eqs. (96) 
and (97). Any higher order radial derivatives ( h uu needs 
only the second radial derivative) can be obtained from 
recurrence relations obtained from the radial Teukolsky 
equation (see [571). 


C. Truncation 


We n ow have all ingredients needed to evaluate Eq. 
(133). However, we need to decide how to truncate the 


infinite sums (over, l, n and to). Since the (retarded) 
metric perturbation is divergent, we do not expect the 
outer sum over / to converge. Instead the “/-modes”, 

1 


/j6± = 

,L UU 


£ •£ c „(m, n, r„) 

i,s (1 54 ) 

M( 2 /w)!r m -4 2 Yi m ( o), 


P±i6) 

ran 2 x l 1 ran ' 


2 


are e xpected to converge to a constant value Sq given by 
(1441. Note that a single Imn mode of h uu requires data 


from limn modes of the Hertz potential with the same 
to and n but different various values of l\. This suggests 
that the most natural way to loop over the modes is to 
first fix values for to and n and calculate all limn modes 
of the Hertz potential up till a certain predetermined 
Zi.max- If we treat the 2 ^^(r 0 )’s with different 

li as forming a vector, then 

1 


M lu i = m»*o ) Cis (m, n ,r 0 ) 




l 


(155) 


defines a linear operator that produces a vector of con¬ 
tributions to the h\i u ’s. The re-expansion of the spin- 
weighted spherical harmonics to ordinary spherical har¬ 
monics, sm A \ 2 , only introduces a coupling between a finite 
number of Zj-modes. The re-expansion for spheroidal to 
spherical harmonics ( 2 on the other hand, causes 
a coupling between an infinite number of Zi-modes. How¬ 
ever, M/i 1 can be seen to decay exponentially away from 
the diagonal. We can thus control the error caused by 
the finite truncation by choosing a maximum value Z max 
for l that is smaller than Zpmax- In practice, the code (at 
each contribution to an mn-mode) determines the largest 
/max such that 


max 

^max 




x + lji 


Mi 


1,1 


Mk _i_2 l \ 

——'J' “ 


(156) 


The final Z max will be the minumum of the /max’s deter¬ 
mined at each mode. We can now sum the contributions 


of all i, s, and j in (154) to obtain all /mn-modes of h u 
with fixed to and n up to Z m ax, 


1 / 


Y,Cis{m,n,r 0 ) 

V i,s,j 

x E^^<^(^o)( 2 6 ron )|r m ^ 2 ) Y lm ( 0). 


(157) 


h,h 


For the final sum over to and n we observe that 


l lmn,± t/— m—n,=h 

't'uu ftj uu 5 


(158) 


and use this to limit our computation to modes with 
n > 0, obtaining the rest by symmetry. Furthermore, to 
calculate AU we only need the orbital average {h l ^ n ). 
For fixed / and to, and large enough n, the magnitude 
of (h l ^ n ’ ± ) decays exponentially with increasing n. This 
leads us to the following truncation procedure. Starting 
at n = 0, we calculate (Zi^ n,± ) for all |m| < Z max . For 
each next n we calculate {h l J£ n,± ), for the same set of to 
as n — 1, except for those to’s for which 


max I 
/</ 

L l max 


(h[ 


mn— 1,± \ 


(£ 0 Lor > 


< e. 


(159) 


In addition, each computation of {h l ^ n,± ) will also pro¬ 
duce an estimate of the relative error on its value. Gen¬ 
erally, this error will increase as we increase n, mainly 
because the integrand (70) becomes more oscillatory. If 
for some n the maximum estimated relative error for 
becomes larger than our target precision e, we 
adjust our target precision upwards to match the esti¬ 
mated error. 


We continue this procedure until (h 1 ^ 71 ^) has converged 
to our target precision e for all to. 


D. Tail estimation 


Adding together all to and n modes, the procedure above 
will give us all (h 1 ^) up to Z max . This allows us to get 
an estimate of ), 


( 16 °) 

1 -0 

However, we expect ( h l uu± ) — (£q ot ) to be order Z -2 . 
Consequently, the missing tail, 

OO 

E (<^>-<£o or >)~^- ( 161 ) 

1-1 +1 tmax 

1 — •'max 1 

With a typical value of Z max = 20 this would gives us only 
two digits of precision. We know the individual /-modes 
to a much higher precision. We can leverage that fact to 
get an estimate of the missing tail. 





















18 


V 

e 

Here 

Akcay et al. [23] 

Barack and Sago [22] 

10 

0.10 

-0.1277540232(10) 

-0.1277540(3) 

-0.1277554(7) 

15 

0.10 

-0.07687063237(5) 

-0.0768706(2) 

-0.0768709(1) 

20 

0.10 

-0.055221659739(6) 

-0.05522166(7) 

-0.05522177(4) 

100 

0.10 

-0.010101234326660(2) 

-0.0101012344(10) 

— 

10 

0.20 

-0.123647888(2) 

-0.123648(3) 

-0.1236493(7) 

15 

0.20 

-0.07431375582(4) 

-0.07431376(9) 

-0.0743140(1) 

20 

0.20 

-0.0534085449572(6) 

-0.05340854(9) 

-0.05340866(4) 

100 

0.20 

-0.0097893279221005(14) 

-0.0097893274(4) 

— 

10 

0.30 

-0.1168019818(2) 

-0.1168020(6) 

-0.1168034(6) 

15 

0.30 

-0.07007684538(10) 

-0.0700768(5) 

-0.0700771(1) 

20 

0.30 

-0.050403774160(6) 

-0.05040377(4) 

-0.05040388(4) 

100 

0.30 

-0.009270280959(2) 

-0.009270281(4) 

— 

10 

0.40 

-0.107220(6) 

-0.107221(2) 

-0.1072221(5) 

15 

0.40 

-0.0641988(4) 

-0.064199(1) 

-0.0641991(1) 

20 

0.40 

-0.0462337(4) 

-0.0462337(9) 

-0.04623383(4) 

100 

0.40 

-0.0085452(4) 

-0.0085453(2) 

— 


Table II. Numerical results for AU from our new code for eccentric orbits around a Schwarzschild a = 0 compared with previous 
frequency domain from Akcay et al. [23] an d time domain results from Barack and Sago }22j . 


Our procedure is as follows. We first calculate the partial 
sums 


st = - ( £ o OT ))- 


1=0 


We can model Sl as an expansion in 1 /L 


(162) 




0° ^-j- 


k =0 


L k ’ 


with the total sum being estimated as 

{h%t) = lim St = S 0 ± . 

L—too 


(163) 


(164) 


We can estimate <Sq by fitting a truncated version of the 


model (163), 


st 


L k 


k =0 


(165) 


to the 5f[. This fit is improved by weighting the data 
point St with their estimated errors. In practice (es¬ 
pecially for low L) this error is dominated by the part 
of (163) that is not included in the truncated model. If 


we assume that all Sj: are order unity, we can estimate 
this error to be approximately L _femax_1 . This estimate 
is good for all but the smallest L, consequently we want 
to exclude those points from the fit. Similarly, the data 
points with the largest L may be so dominated by the 
truncation errors that they would introduce a systematic 
bias in the fit. We thus fit the Sj J data in a window 


To find the optimal values for L min , L max , and fc max , we 
perform the fit for a wide range of (reasonable) choices 
for their values. For each fit we calculate the estimated 


error on . From these fits we choose the ones with the 
lowest estimated error on S*. We then estimate to 
be the median prediction from this set, and estimate the 
error in St to be the maximum of the error estimated 
from the fit and the median deviation of the estimates of 
S* among the best fits. 

This allows us to get a much more accurate estimate of 
(ftjjp). In particular, this estimate can be more accurate 
than the estimated errors on the l -modes with the largest 
l. 


V. RESULTS 


A. Comparison with existing results 


An important verification of the numerical procedure dis¬ 
cussed in this paper is that it reproduces existing results 
from the literature in the relevant limits. The general¬ 
ized redshift invariant for eccentric orbits around a non¬ 
spinning Schwarzschild black hole was first calculated by 
Barack and Sago [22.! using a time domain based code. 
Later, Akcay et al. [23] published more accurate results 
using their frequency domain code. 

In table [TlJ we compare results from our new code with 
these previous publications. For low eccentricities our 
results are consistently within the error bars from |23j 
providing more significant digits. In fact, our more accu¬ 
rate results suggest that the error bars given in [23] are 
slightly pessimistic. Whereas the results from [23] and 
[221 were marginally compatible due to overlapping error 
bars, the tighter error bars on the new data are incon¬ 
sistent with the last digit of [22]. It appears the time 
domain results consistently overestimate the magnitude 
of AU. 













19 



Figure 2. Convergence of the /-modes with increasing n- 
modes (for an orbit with spin a = 0.9, semilatus rectum 
p = 3.32, and eccentricity e = 0.2). Each line represents 
the n-mode contributions to one /-mode with barkers shades 
representing higher values of /. The two horizontal lines rep¬ 
resent the target accuracy and the regularization parameter 
B to which the sum of the n-modes will (approximately) con¬ 
verge for each /-mode. 



Figure 3. The maximal contribution to the /-modes at each 
m and n n-modes (for an orbit with spin a = 0.9, semilatus 
rectum p = 3.32, and eccentricity e = 0.2). The horizontal 
plane represents the value of the accuracy goal. 


The results for [23] were generated using the same com¬ 
puter cluster in Southampton. This gives us the opportu¬ 
nity to compare the computational efficiency of our new 
semi-analytical Mathematica based code, with the older 
C-based frequency domain code which numerically solved 
the linearized field equations. Running on two 16 proces¬ 
sor nodes, the old code took 3915 seconds to calculate all 
the modes for the orbit with p = 10 and e = 0.1. Run¬ 
ning on just one node our code took just 1009 seconds to 
reproduce their result at a similar precision. 


In [23], a procedure for extracting the coefficients of an 
expansion of A U in powers of e and p from their data was 
described, comparing the result to the post-Newtonian 
(PN) predictions. It is relatively easy to produce a similar 
dataset with our new code, but at a higher accuracy, and 


follow their procedure. The results are listed in table III 


We were able to extract all known PN coefficients up to 
e 6 accurate to at least five digits. 


We also obtain fairly accurate estimates for the e 2 and 
e 4 4PN coefficients, which will also have logp terms. The 
e 2 coefficients probably have a similar expansion in tran¬ 
scendental numbers as the 4PN terms of the expansion 
of AU for circular orbits. This means that the coeffi¬ 
cient of e 2 p~ 5 probably is some rational combination of 
log 2, 7 r 2 , and the Euler-Mascheroni constant 7 . With the 
limit number of digits available we are unable to guess 
what this combination is. The coefficient of e 2 p~ 5 logp, 
on other hand is probably a rational number. Based on 
our numerical fit, it is most likely 99/5. We finally get 
some very rough estimates for the coefficients of the e 6 
4PN and e 2 5PN terms, which at least get the sign and 
order of magnitude correctly. 


Shah et al. [31!] previously calculated the redshift invari¬ 


ant for a particle in a circular equatorial orbit around 
a Kerr black hole using similar techniques as here, but 
numerically solving the Teukolksy equation. In Table |IV| 
we compare their results with the results from our new 
code. The results agree to the expected accuracy. In gen¬ 
eral the new code matches the accuracy of the code from 
[301 with exception of some of the strong field modes. 
The accuracy is currently limited by the highest value 
of ^i,max = 30 that our MST Teukolsky solver m can 
currently handle, whereas [30] goes up to /i imax = 80. 
There is no principle obstruction to raising the / lmax 
limit of our code, it simply requires some time and effort 
to further populate some lookup tables used to seed the 
numerical root Ending routines. 


B. Generalized redshift for eccentric orbits in Kerr 


Agreement with existing results in the relevant limits in¬ 
stils some confidence that our code is operating correctly. 
To further check our code we perform some consistency 
checks. 

We first check that the truncated sums have indeed con¬ 
verged to the expected level. Fig. [2] shows the conver¬ 
gence of the sum over n-modes for a sample orbit with 
spin a = 0.9, semilatus rectum p = 3.32, and eccentricity 
e = 0.2. The low /-modes converge much faster than the 
high /-modes, but all /-modes eventually settle into an 
exponential convergence for large n. 

The characteristic two “hump” structure in Fig. [2] is 
the result of a qualitative difference in the convergence 
behaviour of the high and low m-modes. This behaviour 
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Term 

Est. Here 

Coefficient 

Est. [23] 

Exact PN[25j 

Ff 

4 ± 6 • 1(T 12 

4.0002(8) 

4 

e 4 

K 

— 2 ± 4 • HT 10 

-2.00(1) 

-2 

e 6 

P 2 

0 ± 4 • 1(T 9 


0 

p3 

7 ± 6•10“ 9 

7.02(3) 

7 

e 4 

l ±4 - 10“ 7 


1 

4 

eff 

p3 

| ±4 - 1(T 6 


5 

2 

e 2 

bf 

-14.3120980(5) 

-14.5(4) ? 

s -14.31209731 

e 4 

? 

83.38298(7) 


s 83.38296351 

e 6 

p 1 ' 

-36.421(3) 


a -36.42197567 

p 5 

-345.37(5) 

-1500 ±500 


10gP?5 

19.733(5) 

250 ± 100 


e 4 

P 5 

737(4) 



log 

-48.6(4) 



e 6 

P 5 

-185(12) 



log P^w 

7(2) 



e 2 

P 6 

-2000 ± 400 



log pfe 

-40 ± 20 




Table III. Numerical estimates of the coefficients of the ex¬ 
pansion of AU in Schwarzschild as a power series in p -1 and 
e. The first column shows numerical estimates obtained by 
fitting a large data set of eccentric orbits in Schwarzschild 
with 10 < p < 1000 and 0 < e < 0.4 obtained using our new 
code. The second column shows the numerical results from 
Akcay at el. [23]. The last column show the exactly known 
values of the coefficients up to 3PN. 


is shown in Fig. [3j where we see a 3D representation 
of the mn -mode contributions to the Z-modes. At low 
n, the Z-mode is dominated by the contributions from 
the modes with small to. However, these contributions 
decay much more quickly than the large to contributions. 
Consequently, at large values of n the roles are reversed 
and the Z-modes are dominated by m = |Z| contributions. 


A second apparent feature in Fig. [2] is that some of the n- 
mode contributions to the (high) Z-modes are much bigger 
(in magnitude) than the expected final value of the Z- 
mode (as represented by the value of the regularization 
parameter B ). This indicates that the sum over n -modes 
is very oscillatory and has a larger degree of cancellation. 


For this reason our convergence checks in section IV C 
relative to B , aiming for a fixed accuracy. 


are 


Having established convergence of the truncated n-sums, 
we turn our attention to the Z-modes. Fig. [4] shows the 
values of the Z-modes (again for the same reference or¬ 
bit with spin a = 0.9, semilatus rectum p = 3.32, and 
eccentricity e = 0.2) both for the outside (+) and inside 
(—) limits of h uu . As expected, the values of the Z-modes 
do not decay for large Z. Instead the Z-modes converge 
to some constant value. In section [III F| we argued that 
the asymptotic behaviour of our radiation gauge should 



Figure 4. The 1-modes Zi(j(f for an orbit with spin a = 0.9, 
semilatus rectum p = 3.32, and eccentricity e = 0.2. The top 
lines are the values of the retard field before subtraction of 
the regularization parameters. The bottom lines give the val¬ 
ues after subtraction of the regularization parameters. They 
follow the diagonal gridlines which represent l~ 2 decay. 



Figure 5. The difference between the partial sums of the l- 
modes of (ZiJ u ) (excluding completion terms) of an orbit with 
spin a = 0.9, semilatus rectum p = 3.32, and eccentricity e = 
0.2 . As L increases the difference decays as approximately 
Z' 4 (diagonal gridlines). 


match the analytically known values for Lorenz gauge 
regularization parameters (B, C, and D). This can be 
checked numerically, by subtracting the analytical val¬ 
ues of B and C from our numerical data. The results are 
also shown in Fig. [3] After subtraction of B, the Z-modes 
decay as Z -2 . This implies that our data is compatible 
with the value of both the B and the C = 0 terms. This 
agreement exists for both the outside (+) and inside (—) 
limits of h uu . Similar agreement is found for other or¬ 
bits. Besides providing a consistency check of our data, 
this is also the first numerical confirmation of the an¬ 
alytical calculation of the regularization parameters for 
orbits in Kerr for eccentric orbits. 

Not only do the Z-modes of the reconstructed metric ob¬ 
tained from the interior and exterior regions have similar 
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r 

a = —0.9 

a = -0.7 

a = -0.5 

a = 0.5 

a = 0.7 

a = 0.9 

4 

— 

— 

— 

— 

-0.3934469(1) 

-0.325705(1) 






-0.3934452987(3) 

-0.325704499(2) 

5 

— 

— 

— 

-0.3135068(1) 

-0.2776815(1) 

-0.25072625(8) 





-0.3135069374 

-0.2776815518(1) 

-0.2507261858(1) 

6 

— 

— 

— 

-0.23426933(1) 

-0.21718074(2) 

-0.20325499(2) 





-0.2342693592(1) 

-0.2171807422 

-0.2032550190 

7 

— 

— 

— 

-0.188583050(1) 

-0.1788453416 

-0.170547666(3) 





-0.1885830414 

-0.1788453294 

-0.1705476660 

8 

— 

— 

-0.204384597(2) 

-0.158300271(2) 

-0.152122681(2) 

-0.146708039(2) 




-0.2043845850 

-0.1583002707 

-0.1521226786 

-0.1467080374 

10 

-0.151451799(1) 

-0.1456824115(2) 

-0.1404058589(3) 

-0.1201650349(3) 

-0.1171389690(2) 

-0.1143966838 


-0.1514517993 

-0.1456824108 

-0.1404058593 

-0.1201650350 

-0.1171389696 

-0.1143966840 

15 

-0.0833104689 

-0.0819468044 

-0.0806554028 

-0.0751903937 

-0.0742790117 

-0.0734230600 


-0.0833104689(1) 

-0.0819468043(7) 

-0.080655402(6) 

-0.07519039(3) 

-0.074279012(3) 

-0.0734230600(6) 

20 

-0.0581474333 

-0.0575935465 

-0.0570620298 

-0.0547222331 

-0.0543144710 

-0.0539257152 


-0.0581474333 

-0.0575935465 

-0.0570620298 

-0.0547222331 

-0.0543144710 

-0.0539257152 

30 

-0.0365055425 

-0.0363353216 

-0.0361700729 

-0.0354163457 

-0.0352797020 

-0.0351476286 


-0.0365055425 

-0.0363353216 

-0.0361700729 

-0.0354163457 

-0.0352797020 

-0.0351476286 

50 

-0.0210263381 

-0.0209844570 

-0.0209434412 

-0.0207511788 

-0.0207152577 

-0.0206801699 


-0.0210263381 

-0.0209844570 

-0.0209434412 

-0.0207511788 

-0.0207152577 

-0.0206801699 

70 

-0.0147844703 

-0.0147673393 

-0.0147504968 

-0.0146705787 

-0.0146554476 

-0.0146405986 


-0.0147844703 

-0.0147673393 

-0.0147504968 

-0.0146705787 

-0.0146554476 

-0.0146405986 

100 

-0.0102349197 

-0.0102281718 

-0.0102215165 

-0.0101896244 

-0.0101835216 

-0.0101775103 


-0.0102349197 

-0.0102281718 

-0.0102215166 

-0.0101896245 

-0.0101835217 

-0.0101775104 


Table IV. Numerical results for AU for circular orbits in Kerr from our new code, compared with previous results from the 
circular orbit code of Shah et al. 30]. In each cell, the top number in each cell is the result from our new code, and the 
bottom number the result from [30]. (Note that these are not the values published in [30], but values from an updated version 
of that code with some errors fixed.) Parenthetical figures indicate the estimated error on the last displayed digit. When no 
parenthetical figure the error is at least an order of magnitude smaller than the last shown digit. 


P 

e 

a = —0.9 

a = 0.0 

a = 0.9 

Piso + 1 

0.00 

-0.1591225238(2) 

-0.2208475274(2) 

-0.40871659(2) 

Piso + 1 

0.10 

-0.150998651(3) 

-0.209375588(9) 

-0.391205(6) 

Piso + 1 

0.20 

-0.142034175(2) 

-0.195877797(5) 

-0.3635027(6) 

Piso + 1 

0.30 

-0.13127905(3) 

-0.179591180(7) 

-0.327109(4) 

Piso + 1 

0.40 

-0.1182790(4) 

-0.160212(3) 

-0.283764(2) 

Piso + 10 

0.00 

-0.062991300291(4) 

-0.072055057429096(4) 

-0.09093033701(3) 

Piso + 10 

0.10 

-0.061183242(1) 

-0.070241089(2) 

-0.089218442(2) 

Piso + 10 

0.20 

-0.058180786(3) 

-0.066951046(3) 

-0.085406440(1) 

Piso + 10 

0.30 

-0.054043290(2) 

-0.06226898(1) 

-0.0796315937(2) 

Piso + 10 

0.40 

-0.048833176(5) 

-0.0562888674(2) 

-0.072058378(4) 

Piso + 100 

0.00 

-0.009395362239424(0) 

-0.009616383265554(0) 

-0.009942932464012(0) 

Piso + 100 

0.10 

-0.00927472002(8) 

-0.00950017309(4) 

-0.0098334067(4) 

Piso + 100 

0.20 

-0.0089653010(1) 

-0.00918951172(9) 

-0.0095209535(1) 

Piso + 100 

0.30 

-0.00846904194(5) 

-0.00868616413(8) 

-0.00900719351(3) 

Piso + 100 

0.40 

-0.007788207(3) 

-0.00799229975(7) 

-0.0082942072(9) 

Piso + 1000 

0.00 

-0.000993413311181(0) 

-0.000996016937701(0) 

-0.000999595649400(0) 

Piso + 1000 

0.10 

-0.000983179941(7) 

-0.00098584074(1) 

-0.000989496817(4) 

Piso + 1000 

0.20 

-0.000953066056(2) 

-0.000955719551(2) 

-0.000959363369(3) 

Piso + 1000 

0.30 

-0.000903093610(8) 

-0.000905672227(9) 

-0.000909211662(8) 

Piso + 1000 

0.40 

-0.0008332886(4) 

-0.000835722608(2) 

-0.000839062961(5) 


Table V. Numerical data for AU for a selection of orbits with eccentricity e varying from 0.0 (circular) to 0.4, and black hole 
spin a = —0.9 (retrograde), a = 0 (Schwarzschild), or a = 0.9 (prograde). The semi-latus rectum p for each orbit has been 
chosen relative to the semi-latus rectum of the last stable orbit piso at that eccentricity and spin. Parenthetical figures indicate 
the estimated error on the last displayed digit, with (0) indicating that the error is at least an order of magnitude smaller than 
the last shown digit. Note that the data used for the location of piso is accurate only to about 5 digits. For compactness, the 
table does not show the actual value of p used for each orbit. This data is available as supplementary data. 
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difference^ ^R,Lor 



CCK 


CCK 


expansion is obtained by studying the effect of a local 
gauge transformation. It is not clear that this procedure 
will commute with the diagram. In particular, it is very 
well possible that this different route would still induce 
a gauge difference between pointwise values of regular 
metric perturbation in the interior and exterior. 

We conclude by giving values for A U for a selection of 
sample orbits around a spinning black hole in table [V] 
On a single 16 processor node of our computing cluster, 
the computation time for these orbits ranges from about 
2 minutes for the circular orbits to just over 5 and a half 
hours for the 0.4 eccentricity orbits. 


( 'i.Ret,CCK± ) 1 S,CCK±'i d ifferen f e i.R,CCK 
V > n v.u ) f ' b u,v 




fllS 


VI. DISCUSSION AND CONCLUSIONS 


Figure 6. Commutative diagram representing (an idealized 
version of) the CCK metric reconstrunction procedure fol¬ 
lowed in this paper. Commutativity of the diagram suggests 
that we should expect the inside anc outside limits of h^ CCK 
to be the same. However, beware the caveat mentioned in the 
text. 


convergence behaviour, in fact they turn out to converge 
to the same value as can be seen in Fig. [5j At first 
sight, this is some what surprising since the exterior and 
interior metric perturbations are obtained in different ra¬ 
diation gauges. Since the gauge used to obtain the field 
in the interior vacuum region is irregular at infinity, we 
have no reason to expect (h^~) to be the same as Ku + ) 
when we do not include the completion terms. 

An intuitive understanding of why we would expect 
agreement can be gained from considering the diagram 
in Fig. [^representing an (idealized) version of the proce¬ 
dure used to obtain the regular metric perturbation. The 
vertical arrows in this diagram all represent linear opera¬ 
tors, while the horizontal arrows are simple subtractions. 
Hence the diagram commutes. Following the left side 
of the diagram (roughly) corresponds to our procedure; 
Applying the CCK procedure to fields with a non-zero 
source leads to a gauge ambiguity between the values of 
the retarded and singular fields obtained from the inte¬ 
rior and exterior vacuum regions. A similar ambiguity 
is thus expected in their difference, h^ CCK . However, 
if one follows the right-hand side of the diagram, then 
it follows from the fact that h K Lor is a vacuum solution 
of the linearized Einstein equation that there is no am¬ 
biguity in applying the CCK procedure. Consequently, 
commutativity of the diagram implies that there is also 
no gauge ambiguity from other route. 

It must however be noted that the above diagram is not 
exactly the procedure that we follow. In particular, we do 
not obtain the singular field by calculating the singular 
part of f/q and applying the CCK procedure. Instead, the 
Lorenz gauge singular field is obtained as an asymptotic 
power series in Z, and the radiation gauge version of this 


In this paper we have described the first ever numeri¬ 
cal computation of the gravitational perturbations pro¬ 
duced by a small mass to orbiting a Kerr black hole on 
an equatorial eccentric orbit to first order in to. This 
was achieved by applying the CCK metric reconstruction 
to a perturbed Weyl scalar calculated using the semi- 
analytical MST formalism. 

As a demonstration of our code’s capabilities we calcu¬ 
lated the (generalized) redshift gauge invariant, U. In the 
regimes that were previously explored in the literature— 
circular equatorial orbits in Kerr and eccentric orbits in 
Schwarzschild—our code matches the previously avail¬ 
able results. 

The results for eccentric orbits around a Kerr black hole 
are completely new. They have not been calculated nu¬ 
merically, nor are there (to our knowledge) any post- 
Newtonian expansions of the generalized redshift avail¬ 
able to compare against. Comparison against the an¬ 
alytically known regularization parameters verifies that 
the metric perturbations calculated by our code have 
the correct singular structure. It also provides the first 
numerical confirmation of the analytical calculation of 
the mode-sum regularization parameters for gravitational 
perturbations produced by eccentric orbits in Kerr. 

The semi-analytical nature of the MST formalism at the 
basis of our numerical computation means that very high 
accuracies can be achieved at a relatively low cost. In 
particular, the method is unaffected by the numerical 
problems caused by the near degeneracy of ingoing and 
outgoing modes at very low frequencies, which fundamen¬ 
tally limited the accuracy of the Lorenz gauge frequency 
domain code of [23 :. 

The main limiting factor for the accuracy of our code is 
the limit on the number of Z-modes that can be calcu¬ 
lated. This can be traced back to our ability to calculate 
the so-called “renormalized angular momentum” param¬ 
eter v that is needed to evaluate the MST series. This 
involves numerically finding the root of an algebraic equa- 
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tion involving continued fractions. This is only possible 
if a good enough initial guess for v is available. For low 
frequency inodes, this guess is provided by an analytic 
low frequency expansion. However, no such expansion is 
available for higher frequencies. To enable calculations in 
the strong field regime, our code relies on interpolation 
of a table of known high frequency modes. At this stage, 
this table has been populated up to l = 30, limiting our 
code to l < 30 modes. However, with not too much effort 
it should be possible to populate this table to any value 
needed lifting this limitation. 

For high eccentricities a secondary limiting factor to our 
accuracy is the number of n-modes needed to achieve the 
required accuracy. This limitation is inherent to the fre¬ 
quency domain approach. It can be overcome in two prin¬ 
ciple ways: more computation time or better efficiency. 
The first is currently still an option with the most de¬ 
manding orbits in table [V] completing in just over 5 hours 
on a single node of our cluster. On the other hand, our 
code has not been fully optimized yet. For example, our 
code has many internal parameters controlling the inter¬ 
nal working precision of various components of the code. 
These are currently set to “safe” high values to ensure 
loss of internal precision is never an issue. Optimization 
of these values could probably reduce runtime by an order 
of magnitude. Once optimal values for the internal preci¬ 
sion have been established further performance could be 
gained by porting the most CPU intensive parts of the 
code to fixed precision C code. 

With relatively little effort, we were able to substan¬ 
tially improve on the numerical estimates of the post- 
Newtonian expansion of the generalized redshift invariant 
presented in [25] , It should be possible to obtain much 
better results by following the strategy set out in j.32] to 
estimate the post-Newtonian coefficients for the redshift 
from circular orbits. There the calculations were per¬ 
formed in the extreme weak field (with radii ~ 10 2O M), 
where it is much easier to distinguish the logp terms in 
the expansion. In the extreme weak field very good es¬ 
timates for the renormalized angular momentum exist, 
allowing the MST method to solve the field equations 
to essentially arbitrary precision, thousands of digits if 
necessary. With some tuning our new code should also 
work in that regime, and it should be possible to get ex¬ 
tremely accurate results for the generalized redshift in 
the extreme weak field, both for Scliwarzscliild and Kerr 
black holes. 

The method presented in this paper can in principle be 
applied to any quantity constructed from the metric and 
its derivatives. In particular, it should be possible to ob¬ 
tain the full self-force correction to the orbital motion 
(which requires covariant derivatives of the metric). The 
first (and most CPU intensive) stage of the calculation to 


obtain the Hertz potential would remain identical. For 
the new qua ntiti es we would need the relevant counter 
parts of eq. (1321 and table [I] However, the general pro¬ 
cedure for obtaining them should be identical to the pro¬ 
cedure set out in Sec. The occurrence of higher 

derivatives of the metric simply requires the indices i and 
s to have a wider range. In addition, quantities involv¬ 
ing derivatives of the metric will be more singular, and 
may therefore require a smoother choice of the extension 
in (127), and final expression will contain a sum over 
several . However, since the bulk of computational 
time is spent obtaining the Hertz potential, performance 
should be very similar to the implementation for h uu . 

In principle, it should also be possible to apply simi¬ 
lar methodology to obtain the gravitational perturba¬ 
tions produced by a mass on a generic—eccentric and 
inclined—orbit around a Kerr black hole. Our solver for 
the Teukolsky equation m can already handle such or¬ 
bits. Moreover, the general techniques used here (ex¬ 
tended homogeneous solutions, mode-by-mode CCK re¬ 
construction) are equally applicable to the generic case. 
Obviously, the results in sec nn E] and |IIIF| would have 
to be generalized to the generic case. However, we see 
no reason to believe that there are any fundamental ob¬ 
structions (other than complexity) to doing so. 


For generic orbits, computational costs are expected to 
become an issue due to the number of modes that need to 
be calculated in the frequency domain. Some savings can 
however be attained by first considering eccentric inclined 
orbits exhibiting a (low order) y’fTresonance [72]. Access 
to the full self-force on r$-resonant orbits should further 
allow verification for approximate results for the effect of 
such resonances using orbit averaged fluxes EH m]. it 
would also provide valuable input for the open question 
whether motion under the conservative self-force in Kerr 
is integrable m- 
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